

default_zone=zeros(NBG,NB,NSS);

 for j=1:NSS %Today's Exogenous Shock
           for t=2:NB %Today's private debt index;                     
                for i=1:NBG %Today's public Debt shock
                      if Def_func(i,t,j)< Def_func(i,t-1,j)  
                          default_zone(i,t,j)=1;
                      end
                 end
           end
 end
 

T=100000;
cut=500;  % Initial condition dependence
N=1;

addpath('auxiliar')
 load seed
 rng(s);
 slength=30; %Total Length of the Sample (used+Cut)
finalcut=17; % Periods Of Market access before we use (cut)


% 1. Simulation
S=zeros(N,T+cut,N);
S_index=zeros(T+cut,N);

inv_dist=invariant(Prob);
cum_dist=cumsum(inv_dist);
load seed
rng(s);
init0=rand(N,1);
%FInd Invariant Distribution
inva = @(r) find(r<cum_dist,1,'first'); % find the 1st index s.t. r<D(i);
% Now this are your results of the random trials
rR = arrayfun(inva,init0);

for i=1:N
rng(i,'twister');
S=markov(Prob,T+cut+1,rR(i),1:length(Prob));
S_index(:,i)=S';
end

%[S_index,]=markov(Prob,T+cut+1,1,1:length(Prob));
ySIM=zeros(T+cut,N);
yNSIM=zeros(T+cut,N);
KAPPASIM=zeros(T+cut,N);
epsilons_default=zeros(T+cut,N);
epsilons_debt=zeros(T+cut,N);
DefDE_SIM=NaN(T+cut,N);

rng(N+1,'twister');
epsilons_default=rand(T+cut,N);
epsilon_reentry=rand(T+cut,N);
rng(N+2,'twister');
epsilons_debt=rand(T+cut,N);
PVSIM=zeros(T+cut,N);
SSSIM=zeros(T+cut,N);
for i=1:T+cut
   for j=1:N
    ySIM(i,j)=yT(S_index(i,j));
    yNSIM(i,j)=yN(S_index(i,j));
    KAPPASIM(i,j)=KAPPAS(S_index(i,j));
    %PVSIM(i,j)=PVD(S_index(i,j));
   end
end


% Core Interpolations

bpSIM=zeros(T+cut,N);
bSIM=zeros(T+cut,N);
BGSIM=zeros(T+cut,N);
BGpSIM=zeros(T+cut,N);
QSIM=zeros(T+cut,N);
DefSIM=zeros(T+cut,N);
DefSIM_DE=zeros(T+cut,N);
BGp_bigSIM=zeros(T+cut,N,NBG);
qSIM=zeros(T+cut,N);
DebtLimSIM=zeros(T+cut,N);
TAUSIM=zeros(T+cut,N);
BGSIM_Defaulted=zeros(T+cut,N);
bSIM_Defaulted=zeros(T+cut,N);

bSIM_Defaultedm1=zeros(T+cut,N);
bSIM_Defaultedm2=zeros(T+cut,N);
bSIM_Defaultedm3=zeros(T+cut,N);
bSIM_Defaultedm4=zeros(T+cut,N);
bSIM_Defaultedm5=zeros(T+cut,N);
BGSIM_Defaultedm1=zeros(T+cut,N);
BGSIM_Defaultedm2=zeros(T+cut,N);
BGSIM_Defaultedm3=zeros(T+cut,N);
BGSIM_Defaultedm4=zeros(T+cut,N);
BGSIM_Defaultedm5=zeros(T+cut,N);


%Indices
BGiSIM=NaN(T+cut+1,N);
BGpiSIM=NaN(T+cut,N);
zoneSIM=NaN(T+cut,N);
cSIM=zeros(T+cut,N);
priceSIM=zeros(T+cut,N);
SpreadSIM=zeros(T+cut,N);
zeroBGI_ind=find(BGgrid>-0.0001,1);
BadStandindSIM=zeros(T+cut,N);
reentrySIM=zeros(T+cut,N);
BGiSIM(1,:)=zeroBGI_ind;
bSIM(1,:)=1;%0.5;

for i=1:T+cut
    for j=1:N
    bad_stand=BadStandindSIM(i,j);
    yindex=S_index(i,j);
    BGindex=BGiSIM(i,j);
    qSIM(i,j)=q(S_index(i,j));
    bgrid_cc_temp=squeeze(bgrid(BGindex,:));
    zoneSIM(i,j)=interp1(bgrid_def_temp,squeeze(Def_func(BGindex,:,yindex)),bSIM(i,j),'linear','extrap');
    def_temp=min(max(interp1(bgrid_cc_temp,squeeze(Def_func(BGindex,:,yindex)),bSIM(i,j),'linear','extrap'),0),1);
    
   if bad_stand==0 || (bad_stand==1 && theta>epsilon_reentry(i,j)) % WE ARE IN GOOD STANDING or we were in bad standing but we reenter
          if (bad_stand==1 && theta>epsilon_reentry(i,j)) % THIS WAS A REENTRY PERIOD
           reentrySIM(i,j)=1;
          end

        if def_temp>epsilons_default(i,j) %WE ENTER A DEFAULT
            def_temp_SP=min(max(interp1(bgrid_cc_temp,squeeze(Def_funcDE(BGindex,:,yindex)),bSIM(i,j),'linear','extrap'),0),1);
               if def_temp_SP>epsilons_default(i,j)
                DefDE_SIM(i,j)=0; % THE DE WOULD HAVE ALSO ENTERED THIS DEFAULT
            else
                 DefDE_SIM(i,j)=1; % THE DE WOULD NOT HAVE ALSO ENTERED THIS DEFAULT
               end

            BadStandindSIM(i+1,j)=1; % TOMORROW WILL BE IN BAD STANDING
            DefSIM(i,j)=true(1);            
            BGpSIM(i,j)=0;
            BGSIM(i,j)=0;
            BGSIM_Defaulted(i,j)=BGgrid(BGindex);
            BGpiSIM(i,j)=zeroBG_ind;
            BGiSIM(i,j)=zeroBG_ind;
            BGindex=zeroBG_ind;
            bgrid_def_temp=squeeze(bgrid(BGindex,:));
            bpSIM(i,j)=interp1(bgrid_def_temp,squeeze(bp_bad(:,yindex)),bSIM(i,j),'linear','extrap');
            cSIM(i,j)=max(interp1(bgrid_def_temp,squeeze(c_big_bad(:,yindex)),bSIM(i,j),'linear','extrap'),1e-8);
            %priceSIM(i,j)=interp1(bgrid_def_temp,squeeze(price(zeroBG_ind,:,yindex,zeroBG_ind)),bSIM(i,j),'linear','extrap');
            priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
            QSIM(i,j)=interp1(bgrid_def_temp,squeeze(Qgov(zeroBG_ind,:,yindex,zeroBG_ind)),bSIM(i,j),'linear','extrap');
            TAUSIM(i,j)=(interp1(bgrid_def_temp,squeeze(Tau_bad(:,yindex)),bSIM(i,j),'linear','extrap'));
            SpreadSIM(i,j)=0;
            DebtLimSIM(i,j)=(KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))/qSIM(i,j);
            bpSIM(i,j)=min(bpSIM(i,j), DebtLimSIM(i,j));
            bSIM_Defaulted(i,j)=bSIM(i,j);
            bSIM_Defaultedm1(i,j)=bSIM(i-1,j);
            bSIM_Defaultedm2(i,j)=bSIM(i-2,j);
            bSIM_Defaultedm3(i,j)=bSIM(i-3,j);
            bSIM_Defaultedm4(i,j)=bSIM(i-4,j);
            bSIM_Defaultedm5(i,j)=bSIM(i-5,j);
            BGSIM_Defaultedm1(i,j)=BGSIM(i-1,j);
            BGSIM_Defaultedm2(i,j)=BGSIM(i-2,j);
            BGSIM_Defaultedm3(i,j)=BGSIM(i-3,j);
            BGSIM_Defaultedm4(i,j)=BGSIM(i-4,j);
            BGSIM_Defaultedm5(i,j)=BGSIM(i-5,j);

        else
            DefSIM(i,j)=false(1);  
            BadStandindSIM(i+1,j)=0; % TOMORROW WILL BE IN GOOD STANDING
            for jj=1:NBG
                gtemp=min(max(interp1(bgrid_cc_temp,squeeze(G_big(BGindex,:,yindex,jj)),bSIM(i,j),'linear','extrap'),0),1);
                BGp_bigSIM(i,j,jj)=gtemp+ BGp_bigSIM(i,j,max(jj-1,1)); %We build the Cumulative Distribution
                if epsilons_debt(i,j)<=BGp_bigSIM(i,j,jj) && isnan(BGpiSIM(i,j))==1 % We find the spot in which the shoc crosses the cumulative
                    BGpiSIM(i,j)=jj;
                end
            end             
            BGpSIM(i,j)=BGgrid(BGpiSIM(i,j));
            BGSIM(i,j)=BGgrid(BGindex);
            bpSIM(i,j)=interp1(bgrid_cc_temp,squeeze(bp(BGindex,:,yindex,BGpiSIM(i))),bSIM(i,j),'linear','extrap');
            cSIM(i,j)=max(interp1(bgrid_cc_temp,squeeze(c_big(BGindex,:,yindex,BGpiSIM(i))),bSIM(i,j),'linear','extrap'),1e-8);
            priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
          
            QSIM(i,j)=interp1(bgrid_cc_temp,squeeze(Qgov(BGindex,:,yindex,BGpiSIM(i))),bSIM(i,j),'linear','extrap');
            TAUSIM(i,j)=(interp1(bgrid_cc_temp,squeeze(Tau1p(BGindex,:,yindex,BGpiSIM(i))),bSIM(i,j),'linear','extrap'));
            SpreadSIM(i,j)=(1+(delta)/QSIM(i,j)-delta)-(1+rbar);
            DebtLimSIM(i,j)=(KAPPASIM(i,j)*((priceSIM(i,j)*yNSIM(i,j))+ySIM(i,j)))/qSIM(i,j);
            bpSIM(i,j)=min(bpSIM(i,j), DebtLimSIM(i,j));
            BGSIM_Defaulted(i,j)=NaN;
        end

   else % WE ARE IN BAD STANDING and we dont enter
          BadStandindSIM(i+1,j)=1; % TOMORROW WILL BE IN BAD STANDING
            DefSIM(i,j)=true(1);            
            BGpSIM(i,j)=0;
            BGSIM(i,j)=0;
            BGSIM_Defaulted(i,j)=NaN;
            BGpiSIM(i,j)=zeroBG_ind;
            BGiSIM(i,j)=zeroBG_ind;
            BGindex=zeroBG_ind;
            bgrid_def_temp=squeeze(bgrid(BGindex,:));
            bpSIM(i,j)=interp1(bgrid_def_temp,squeeze(bp_bad(:,yindex)),bSIM(i,j),'linear','extrap');
            cSIM(i,j)=max(interp1(bgrid_def_temp,squeeze(c_big_bad(:,yindex)),bSIM(i,j),'linear','extrap'),1e-8);
          
            priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
            QSIM(i,j)=interp1(bgrid_def_temp,squeeze(Qgov(zeroBG_ind,:,yindex,zeroBG_ind)),bSIM(i,j),'linear','extrap');
            TAUSIM(i,j)=(interp1(bgrid_def_temp,squeeze(Tau_bad(:,yindex)),bSIM(i,j),'linear','extrap'));
            SpreadSIM(i,j)=0;
            DebtLimSIM(i,j)=(KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))/qSIM(i,j);
            bpSIM(i,j)=min(bpSIM(i,j), DebtLimSIM(i,j));

   end
    BGiSIM(i+1,j)=BGpiSIM(i,j);
    bSIM(i+1,j)=bpSIM(i,j);
    end
end

% Computation

bSIM=bSIM(1:end-1,:);%

% Current Accpunt Calculation
CA_SIM=bpSIM(2:end,:)-bpSIM(1:end-1,:);  % Current Account of the Private Sector

cSIM=cSIM(cut+1:T+cut,:);
ySIM=ySIM(cut+1:T+cut,:);
yNSIM=yNSIM(cut+1:T+cut,:);
S_index=S_index(cut+1:T+cut,:);
DefSIM=DefSIM(cut+1:T+cut,:); 
BGiSIM=BGiSIM(cut+1:T+cut,:);
BGpiSIM=BGpiSIM(cut+1:T+cut,:); 
BGpSIM=BGpSIM(cut+1:T+cut,:); 
bpSIM=bpSIM(cut+1:T+cut,:);
BGSIM=BGSIM(cut+1:T+cut,:);
bSIM=bSIM(cut+1:T+cut,:);
QSIM=QSIM(cut+1:T+cut,:);
KAPPASIM=KAPPASIM(cut+1:T+cut,:);
BGp_bigSIM=BGp_bigSIM(cut+1:T+cut,:,:);
epsilons_debt=epsilons_debt(cut+1:T+cut,:);
epsilons_default=epsilons_default(cut+1:T+cut,:);
DebtLimSIM=DebtLimSIM(cut+1:T+cut,:);
qSIM=qSIM(cut+1:T+cut,:);
BGSIM_Defaulted=BGSIM_Defaulted(cut+1:T+cut,:);
CA_SIM=CA_SIM(cut:T+cut-1,:);
SpreadSIM=SpreadSIM(cut:T+cut-1,:);
BadStandindSIM=BadStandindSIM(cut:T+cut-1,:);
priceSIM=priceSIM(cut+1:T+cut,:);   
TAUSIM=TAUSIM(cut+1:T+cut,:);
reentrySIM=reentrySIM(cut+1:T+cut,:);
bSIM_Defaulted=bSIM_Defaulted(cut+1:T+cut,:);
zoneSIM=zoneSIM(cut+1:T+cut,:);
DefDE_SIM=DefDE_SIM(cut+1:T+cut,:);
     
            bSIM_Defaultedm1=bSIM_Defaultedm1(cut+1:T+cut,:);
            
            bSIM_Defaultedm2=bSIM_Defaultedm2(cut+1:T+cut,:);
            bSIM_Defaultedm3=bSIM_Defaultedm3(cut+1:T+cut,:);
            bSIM_Defaultedm4=bSIM_Defaultedm4(cut+1:T+cut,:);
            bSIM_Defaultedm5=bSIM_Defaultedm5(cut+1:T+cut,:);
            
            BGSIM_Defaultedm1=BGSIM_Defaultedm1(cut+1:T+cut,:);
            BGSIM_Defaultedm2=BGSIM_Defaultedm2(cut+1:T+cut,:);
            BGSIM_Defaultedm3=BGSIM_Defaultedm3(cut+1:T+cut,:);
            BGSIM_Defaultedm4=BGSIM_Defaultedm4(cut+1:T+cut,:);
            BGSIM_Defaultedm5=BGSIM_Defaultedm5(cut+1:T+cut,:);


VSIM=zeros(T,N);
bindSIM=zeros(T,N);
TransSIM=zeros(T,N);
TightSIM=zeros(T,N);
NUSSIM=zeros(T,N);
spreadSIM=zeros(T,N);
priceSIM=zeros(T,N);

for i=1:T
    for j=1:N
    bgrid_temp=squeeze(bgrid(BGiSIM(i,j),:)); 
    VSIM(i,j)=interp1(bgrid_temp,squeeze(V(BGiSIM(i,j),:,S_index(i))),bSIM(i,j),'linear','extrap');

 
   
    if DefSIM(i,j)==0
        TransSIM(i,j)=QSIM(i,j)*(BGpSIM(i,j)-(1-delta)*BGSIM(i,j))-delta*BGSIM(i,j);
        spreadSIM(i,j)=(delta-delta*QSIM(i,j))/QSIM(i,j)-rbar;
        
    else
        TransSIM(i,j)=0;
    end
    priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
    DebtLimSIM(i,j)=(KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))/qSIM(i,j);
   if abs(bpSIM(i,j)-(1/qSIM(i,j))*KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))<1e-3
            bindSIM(i,j)=1;
            TightSIM(i,j)=bpSIM(i,j)-(1/qSIM(i,j))*KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j));          
   end

    end
end

BGSSm3_SIM=NaN(T,N);
BGSSm2_SIM=NaN(T,N);
BGSSm1_SIM=NaN(T,N);
BGSS_SIM=NaN(T,N);
BGSSp1_SIM=NaN(T,N);
BGSSp2_SIM=NaN(T,N);
BGSSp3_SIM=NaN(T,N);


bSSm3_SIM=NaN(T,N);
bSSm2_SIM=NaN(T,N);
bSSm1_SIM=NaN(T,N);
bSS_SIM=NaN(T,N);
bSSp1_SIM=NaN(T,N);
bSSp2_SIM=NaN(T,N);
bSSp3_SIM=NaN(T,N);


Y_SIM=priceSIM.*yNSIM+ySIM;


% SUDDEN STOPS
CAo_SIM=CA_SIM./Y_SIM;  % Current Account
% Current Account and SS
SS=zeros(T,1);
SS(CAo_SIM<(mean(CAo_SIM)-2*std(CAo_SIM)))=1; 

for i=1:T-6
    for j=1:N
        if SS(i+2,j)==1 
        BGSSm3_SIM(i,j)=BGSIM(i,j);
        BGSSm2_SIM(i,j)=BGSIM(i+1,j);
        BGSSm1_SIM(i,j)=BGSIM(i+2,j);
        BGSS_SIM(i,j)=BGSIM(i+3,j);
        BGSSp1_SIM(i,j)=BGSIM(i+4,j);
        BGSSp2_SIM(i,j)=BGSIM(i+5,j);
        BGSSp3_SIM(i,j)=BGSIM(i+6,j);

        bSSm3_SIM(i,j)=bSIM(i,j);
        bSSm2_SIM(i,j)=bSIM(i+1,j);
        bSSm1_SIM(i,j)=bSIM(i+2,j);
        bSS_SIM(i,j)=bSIM(i+3,j);
        bSSp1_SIM(i,j)=bSIM(i+4,j);
        bSSp2_SIM(i,j)=bSIM(i+5,j);
        bSSp3_SIM(i,j)=bSIM(i+6,j);
        end
    end
end


CSIM=(omega*cSIM.^(-ita)+(1-omega)*yNSIM.^(-ita)).^(-1/ita);
meanY=mean(priceSIM.*yNSIM+ySIM);

 mean_default=mean(DefSIM);
cte_debtfactor=delta/(1-(1-delta)/R) ;
amean_PubDebt=cte_debtfactor*mean(BGpSIM)/mean(priceSIM.*yNSIM+ySIM);
amean_PvdDebt=mean(bpSIM)/mean(priceSIM.*yNSIM+ySIM);
% Private Crisis
BQSIM=(BGpSIM-(1-delta)*BGSIM).*QSIM;
avSpr=mean(spreadSIM); 

PvDebt_to_GDP=bpSIM./(priceSIM.*yNSIM+ySIM);
PubDebt_to_GDP=cte_debtfactor*BGpSIM./(priceSIM.*yNSIM+ySIM);
avPvD=mean(PvDebt_to_GDP);
avPub=mean(PubDebt_to_GDP);
 NFA=PubDebt_to_GDP+PvDebt_to_GDP;
      BGSIM_Defaulted_SP=cte_debtfactor*BGSIM_Defaulted(~isnan(BGSIM_Defaulted))/meanY;
                  bSIM_Defaulted_SP=bSIM_Defaulted(~isnan(BGSIM_Defaulted))/meanY;
bSIM_Defaultedm1_SP=bSIM_Defaultedm1(~isnan(BGSIM_Defaulted))/meanY;
bSIM_Defaultedm2_SP=bSIM_Defaultedm2(~isnan(BGSIM_Defaulted))/meanY;
bSIM_Defaultedm3_SP=bSIM_Defaultedm3(~isnan(BGSIM_Defaulted))/meanY;
bSIM_Defaultedm4_SP=bSIM_Defaultedm4(~isnan(BGSIM_Defaulted))/meanY;
bSIM_Defaultedm5_SP=bSIM_Defaultedm5(~isnan(BGSIM_Defaulted))/meanY;


BGSIM_Defaultedm1_SP=cte_debtfactor*BGSIM_Defaultedm1(~isnan(BGSIM_Defaulted))/meanY;
BGSIM_Defaultedm2_SP=cte_debtfactor*BGSIM_Defaultedm2(~isnan(BGSIM_Defaulted))/meanY;
BGSIM_Defaultedm3_SP=cte_debtfactor*BGSIM_Defaultedm3(~isnan(BGSIM_Defaulted))/meanY;
BGSIM_Defaultedm4_SP=cte_debtfactor*BGSIM_Defaultedm4(~isnan(BGSIM_Defaulted))/meanY;
BGSIM_Defaultedm5_SP=cte_debtfactor*BGSIM_Defaultedm5(~isnan(BGSIM_Defaulted))/meanY;


BGSSm3_SP=cte_debtfactor*BGSSm3_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSSm2_SP=cte_debtfactor*BGSSm2_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSSm1_SP=cte_debtfactor*BGSSm1_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSS_SP=cte_debtfactor*BGSS_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSSp1_SP=cte_debtfactor*BGSSp1_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSSp2_SP=cte_debtfactor*BGSSp2_SIM(~isnan(BGSSm3_SIM))/meanY;
BGSSp3_SP=cte_debtfactor*BGSSp3_SIM(~isnan(BGSSm3_SIM))/meanY;


bSSm3_SP=bSSm3_SIM(~isnan(bSSm3_SIM))/meanY;
bSSm2_SP=bSSm2_SIM(~isnan(bSSm3_SIM))/meanY;
bSSm1_SP=bSSm1_SIM(~isnan(bSSm3_SIM))/meanY;
bSS_SP=bSS_SIM(~isnan(bSSm3_SIM))/meanY;
bSSp1_SP=bSSp1_SIM(~isnan(bSSm3_SIM))/meanY;
bSSp2_SP=bSSp2_SIM(~isnan(bSSm3_SIM))/meanY;
bSSp3_SP=bSSp3_SIM(~isnan(bSSm3_SIM))/meanY;

DefDE_SP=DefDE_SIM(~isnan(DefDE_SIM));

fprintf('The long-run probability of crisis in SP is: ');
fprintf('%5.4f \n ',  mean(SS)*100);


fprintf('The long-run probability of being in the shrinking default zone in SP is: ');
fprintf('%5.4f \n ',  mean(zoneSIM)*100);

fprintf('Total number of defaults in SP is: ');
fprintf('%5.4f \n ',  size(DefDE_SP));

fprintf('Total number of defaults in SP that DE would have avoided is: ');
fprintf('%5.4f \n ',  sum(DefDE_SP));



